home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 21 / Mac Magazin and MacEasy Magazine CD - Issue 21.iso / Wissenschaft & Technik / yorick12vr1-nofpu folder / include / bowtie.i < prev    next >
Text File  |  1995-07-26  |  4KB  |  139 lines
  1. /*
  2.    BOWTIE.I
  3.    Detect bowties, boomerangs, and other topological oddities in
  4.    a quadrilateral mesh.
  5.  
  6.    $Id$
  7.  */
  8. /*    Copyright (c) 1994.  The Regents of the University of California.
  9.                     All rights reserved.  */
  10.  
  11. func bowtie(rt, zt, ireg)
  12. /* DOCUMENT map= bowtie(rt, zt)
  13.          or map= bowtie(rt, zt, ireg)
  14.      returns a "bowtie map" for the quadrilateral mesh defined by
  15.      RT, ZT, and (optionally) IREG.  If IREG is present, it should be
  16.      an integer array of the same dimensions as RT and ZT; its first
  17.      row and column are ignored, otherwise each non-zero element of
  18.      IREG marks an existing zone in the mesh.  (An IREG with one fewer
  19.      row and column than RT and ZT will also be accepted.)  If IREG
  20.      is omitted, every zone is presumed to exist.
  21.  
  22.      The returned MAP is a 2-D integer array with one fewer row and
  23.      column than RT and ZT.  It's values have the following meanings:
  24.  
  25.           2   marks a convex zone with positive area
  26.       1   marks a concave (boomerang) zone with positive area
  27.       0   marks a bowtied zone
  28.      -1   marks a concave (boomerang) zone with negative area
  29.      -2   marks a convex zone with negative area
  30.      -9   marks a non-existent zone
  31.  
  32.      Use the nbow function to print the results.
  33.  
  34.    SEE ALSO: nbow
  35.  */
  36. {
  37.   dims= dimsof(rt);
  38.   d= dimsof(zt);
  39.   if (d(1)!=dims(1) || dims(1)!=2 || anyof(dims!=d) ||
  40.       dims(2)<2 || dims(3)<2)
  41.     error, "rt, zt are not quadrilateral mesh coordinates";
  42.  
  43.   dz1= zt(dif,1:-1);    dr1= rt(dif,1:-1);
  44.   dz2= zt(2:0,dif);     dr2= rt(2:0,dif);
  45.   dz3= zt(dif,2:0);     dr3= rt(dif,2:0);
  46.   dz4= zt(1:-1,dif);    dr4= rt(1:-1,dif);
  47.  
  48.   a= array(0.0, 4, dimsof(dz1));
  49.   a(1,,)= dz1*dr4 - dz4*dr1;
  50.   a(2,,)= dz1*dr2 - dz2*dr1;
  51.   a(3,,)= dz3*dr2 - dz2*dr3;
  52.   a(4,,)= dz3*dr4 - dz4*dr3;
  53.  
  54.   map= long(sign(a))(sum,,)/2;
  55.  
  56.   if (is_void(ireg)) return map;
  57.  
  58.   if ((d=dimsof(ireg))(1)==2 && d(2)==dims(2) && d(3)==dims(3)) {
  59.     ireg= ireg(2:0, 2:0);
  60.   } else if (d(1)!=2 || d(2)!=dims(2)-1 || d(3)!=dims(3)-1) {
  61.     error, "ireg not conformable with rt, zt arrays";
  62.   }
  63.  
  64.   if (allof(ireg)) return map;
  65.  
  66.   map(where(!ireg))= -9;
  67.   return map;
  68. }
  69.  
  70. local nbow_negative;
  71. /* DOCUMENT nbow_negative
  72.      default value is 0
  73.      set to 1 to reverse the sense of positive area for the nbow function
  74.  */
  75. nbow_negative= 0;
  76.  
  77. func nbow(rt, zt, ireg, all=)
  78. /* DOCUMENT nbow, map
  79.          or nbow, file
  80.          or nbow, rt, zt
  81.      or nbow, rt, zt, ireg
  82.      prints information about topological oddities in a mesh.
  83.      MAP is a bowtie map as returned by the bowtie function.
  84.      FILE is a binary file containing rt, zt, and ireg arrays.
  85.      RT, ZT and IREG are 2-D arrays defining a quadrilateral mesh.
  86.  
  87.      The information printed includes the zone index (corner with
  88.      the largest indices) of zones which are concave (boomerangs)
  89.      or bowtied, and of zones with negative area.  You can set
  90.      the global variable nbow_negative to 1 to reverse the default
  91.      sense of positive area.  By default, only the first 10 zones
  92.      in each category are printed; use the all=1 keyword argument
  93.      to print a complete (and maybe very long) list.
  94.  
  95.    SEE ALSO: bowtie
  96.  */
  97. {
  98.   if (!is_void(zt)) map= bowtie(rt, zt, ireg);
  99.   else if (is_stream(rt)) {
  100.     f= rt;
  101.     vars= get_vars(f);
  102.     if ((vars(1) && anyof(*vars(1) == "ireg")) ||
  103.     (vars(2) && anyof(*vars(2) == "ireg"))) ireg= f.ireg;
  104.     else ireg= [];
  105.     map= bowtie(f.rt, f.zt, ireg);
  106.   } else {
  107.     map= rt;
  108.   }
  109.  
  110.   if (nbow_negative) map= -map;
  111.  
  112.   list= where2(abs(map)<=2 & map<0);
  113.   _nbow_print, list, "Found %ld zones with negative areas:\n";
  114.  
  115.   list= where2(map==0);
  116.   _nbow_print, list, "Found %ld bowtied zones:\n";
  117.  
  118.   list= where2(map==1);
  119.   _nbow_print, list, "Found %ld concave (boomerang) zones:\n";
  120.  
  121.   list= where(map==2);
  122.   write, format="Found %ld convex, positive area zones\n", numberof(list);
  123. }
  124.  
  125. func _nbow_print(list, format)
  126. {
  127.   n= numberof(list);
  128.   if (n) {
  129.     list++;
  130.     n/= 2;
  131.     write, format=format, n;
  132.     if (!all && n>10) np= 10;
  133.     else np= n;
  134.     write, format="     k= %3ld   l= %3ld\n", list(1,1:np), list(2,1:np);
  135.     if (np<n)
  136.       write, format="     (%ld more, use all=1 keyword to see them)\n", n-np;
  137.   }
  138. }
  139.